In [1]:
from util import *
from glob import glob

/usr/local/lib/python3.10/dist-packages/geopandas/_compat.py:123: UserWarning: The Shapely GEOS version (3.11.1-CAPI-1.17.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(
/home/nyou045/git/CliffErosion_Projections/util.py:3: UserWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas will still use PyGEOS by default for now. To force to use and test Shapely 2.0, you have to set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas:

import os
os.environ['USE_PYGEOS'] = '0'
import geopandas

In a future release, GeoPandas will switch to using Shapely by default. If you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html).
  import geopandas as gpd  # Geospatial data
In [2]:
df = load_AOIs()
df
Out[2]:
Taranaki AOI SSP 4.5 (p50) SSP 4.5 (p83) SSP 8.5 (p50) SSP 8.5 (p83) Rate SSP 4.5 (p50) Rate SSP 4.5 (p83) Rate SSP 8.5 (p50) Rate SSP 8.5 (p83) match match_score
7 NORTH TongaporutuRiver 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 TongapurutuRiverCliffs 93.750000
11 SOUTH HangatahuaRiver_South 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 HangatahuRiver_South 97.435897
21 SOUTH Rahotu 0.58 0.78 0.84 1.10 0.0058 0.0078 0.0084 0.0110 Rahotu 100.000000
20 SOUTH Pihama 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Pihama 100.000000
19 SOUTH OpunakeBeach 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OpunakeBeachCliffs 100.000000
18 SOUTH OhaweBeach 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 OhaweBeach 100.000000
17 SOUTH Oeo 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Oeo 100.000000
16 SOUTH Manutahi 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Manutahi 100.000000
15 SOUTH ManaBay 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 ManaBayCliffs 100.000000
14 SOUTH KaupokonuiBeach 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 KaupokonuiBeach 100.000000
13 SOUTH Kakaramea 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Kakaramea 100.000000
12 SOUTH Hawera_WaihiBeach 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Hawera_WaihiBeach 100.000000
0 NORTH Mohakatino_River 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 MohakatinoRiver 100.000000
10 SOUTH CapeEgmont 0.58 0.78 0.84 1.10 0.0058 0.0078 0.0084 0.0110 CapeEgmont 100.000000
9 NORTH Waitara 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Waitara 100.000000
8 NORTH UrenuiRiver_North 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 UrenuiRiverNorthCliffs 100.000000
6 NORTH PariokariwaPoint 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 PariokariwaPointCliffs 100.000000
5 NORTH Onaero 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OnaeroCliff 100.000000
4 NORTH Oakura_South 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OakuraSouth 100.000000
3 NORTH Oakura 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Oakura 100.000000
2 NORTH NewPlymouth_Airport 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 NewPlymouthAirport 100.000000
1 NORTH NewPlymouth 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 NewPlymouthCliffs 100.000000
22 SOUTH WainuiBeach 0.57 0.78 0.83 1.09 0.0057 0.0078 0.0083 0.0109 WainuiBeach 100.000000
23 SOUTH WaipipiBeach 0.57 0.78 0.83 1.09 0.0057 0.0078 0.0083 0.0109 WaipipiBeachCliffs 100.000000
In [3]:
site = df.match.sample(1).iloc[0]
site
Out[3]:
'UrenuiRiverNorthCliffs'
In [4]:
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
if site == "ManaBayCliffs":
  print("Flipping")
  for k, v in transect_metadata.items():
    transect_metadata[k]["Azimuth"] = v["Azimuth"] + 180
In [5]:
results = predict(gdf, transect_metadata)
results
Out[5]:
TransectID BaselineID group Year ocean_point linear_model_point sqrt_model_point BH_model_point Sunamura_model_point
0 5 1 1 2100 POINT (1720633.1100425983 5684136.280603779) POINT (1720854.3339775638 5683672.05802439) POINT (1720868.1061666948 5683643.158072117) POINT (1720869.6189754275 5683639.983551366) POINT (1720854.3339775638 5683672.05802439)
1 6 1 1 2100 POINT (1720642.0801050563 5684139.860357197) POINT (1720862.6898028378 5683676.926710536) POINT (1720873.4603526713 5683654.325483295) POINT (1720877.9748007017 5683644.852237511) POINT (1720862.6898028378 5683676.926710536)
2 7 1 1 2100 POINT (1720650.575075268 5684144.8879364645) POINT (1720871.3119035014 5683681.687515657) POINT (1720881.2773631075 5683660.775712584) POINT (1720886.5969013653 5683649.613042634) POINT (1720871.3119035014 5683681.687515657)
3 8 1 1 2100 POINT (1720659.1581115604 5684150.087639844) POINT (1720879.2862932794 5683688.164420146) POINT (1720888.6494842013 5683668.51643476) POINT (1720894.571291143 5683656.089947121) POINT (1720879.2862932794 5683688.164420146)
4 9 1 1 2100 POINT (1720667.0028565435 5684155.18472246) POINT (1720886.9824249977 5683693.57335705) POINT (1720897.615215984 5683671.261206882) POINT (1720902.2674228614 5683661.498884027) POINT (1720886.9824249977 5683693.57335705)
... ... ... ... ... ... ... ... ... ...
314 414 25 31 2100 POINT (1726658.0317596905 5689964.909653514) POINT (1727097.8887274514 5689705.400302266) POINT (1727116.2523344012 5689694.5660355035) POINT (1727128.490082377 5689687.345938357) POINT (1727097.8887274514 5689705.400302266)
315 415 25 31 2100 POINT (1726662.3798315164 5689973.425984194) POINT (1727101.3898836626 5689714.932497662) POINT (1727116.3902041595 5689706.100159711) POINT (1727132.0069823388 5689696.904845344) POINT (1727101.3898836626 5689714.932497662)
316 416 25 31 2100 POINT (1726666.1429697142 5689982.1815043455) POINT (1727103.4452041262 5689725.309927042) POINT (1727114.616595727 5689718.747845362) POINT (1727134.0811645247 5689707.314346753) POINT (1727103.4452041262 5689725.309927042)
317 417 25 31 2100 POINT (1726669.315968925 5689991.279733421) POINT (1727107.1794726513 5689734.6948276665) POINT (1727119.1615009448 5689727.67344412) POINT (1727137.834261176 5689716.731339142) POINT (1727107.1794726513 5689734.6948276665)
318 418 25 31 2100 POINT (1726672.9711513557 5690000.091062417) POINT (1727113.035012376 5689742.835461913) POINT (1727129.9653122795 5689732.938229137) POINT (1727143.70859541 5689724.904084851) POINT (1727113.035012376 5689742.835461913)

319 rows × 9 columns

In [6]:
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
poly = prediction_results_to_polygon(results)
In [7]:
poly.explore(tiles="Esri.WorldImagery")
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [8]:
site = df.match.sample(1).iloc[0]
site = "CapeEgmont"
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
m = gpd.read_file(f"Shapefiles/{site}_TransectLines.shp").explore("TransectID", tiles="Esri.WorldImagery", max_zoom=22)
results = predict(gdf, transect_metadata)
results["geometry"] = results["linear_model_point"]
pd.concat([gdf, results])[["Year", "geometry"]].explore("Year", m=m)
gpd.GeoSeries(results.ocean_point, crs=2193).explore(m=m, color="blue")
#gpd.GeoSeries(poly, crs=2193).explore(color="pink", m=m)
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
polygon = prediction_results_to_polygon(results)
polygon.explore(m=m)
m

CapeEgmont
Out[8]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [9]:
run_all_parallel()

  0%|          | 0/24 [00:00<?, ?it/s]
Flipping
In [10]:
def read_file(f):
  df = gpd.read_file(f)
  df["filename"] = f
  return df
samples = pd.concat(read_file(f) for f in glob("Projected_Shoreline_Polygons/*.shp"))
samples.explore("filename", tiles="Esri.WorldImagery", style_kwds={"fill": False})
Out[10]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [11]:
# Compare different SLR rates
site = df.match.sample(1).iloc[0]
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
low_proj_slr = predict(gdf, transect_metadata, Proj_SLR=0.007)
high_proj_slr = predict(gdf, transect_metadata, Proj_SLR=0.015)
m = gpd.GeoSeries(low_proj_slr.sqrt_model_point, crs=2193).explore(color="blue", tiles="Esri.WorldImagery")
gpd.GeoSeries(high_proj_slr.sqrt_model_point, crs=2193).explore(color="green", m=m)

OpunakeBeachCliffs
Out[11]:
Make this Notebook Trusted to load map: File -> Trust Notebook